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Method and App aratus for Evaluating Polynomials and Rational Functions 

Abstract 

Disclosed herein are a computer-processing method and apparatus for computing values of 
polynomials or rational functions. A mathematical software library can advantageously embody the 
concepts of this invention. The method can be adapted to compute values for 
non-elementary special functions, for example ERF, ERFC, LGAMMA, and Bessel functions. The 
steps for polynomial evaluation include presenting input data that includes coefficients of 
polynomial p(x) 9 x, a predetermined x„ and p(x$, building polynomial c(x) having coefficients so 
that polynomial p(x) is expressible as: p(x) = p(xi) + -c(x), determining each coefficient of 

polynomial c(x), determining a value of polynomial c(x) % and constructing a value of polynomial 
p(x) by determining: p(x) =/7(*/) + {*-*/} c(x). The method can be adapted for providing a value 
for a rational function r(x) = p(x)/q(x), which is a ratio of a numerator polynomial p(x) and a 
denominator polynomial q(x). 
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Method and Apparatus for Evaluating Polynomials and Rational Functions 
Field of the Invention 

The present invention relates to a method and apparatus for computing values of 
mathematical expressions, and more specifically to a computer processing method and computer 
5 apparatus for computing binary representations of numerical values of polynomials and rational 
functions. 

Background of the Present Invention 

Polynomials and rational functions, which are quotients of polynomials, are used, for 

example, by various branches of applied science for determining numerical values of mathematical 

io expressions for modelling a property of a physical system, for example a rate of air flow over an 

airfoil surface. Polynomials can be used to approximate complicated mathematical expressions. 

Various methods, for example Horner's Rule that was disclosed in 1819, are used for computing 

values of polynomials, and provide a degree of accuracy that may not be acceptable in certain 

situations, depending on the particular polynomial and the accuracy required. G.E. Forsythe et al in 

15 "Computer Methods for Mathematical Computations" Prentice-Hall (1977) discloses, in Section 

4.2, using Horner's Rule for computing a numerical value of a polynomial. 

Steps for computing binary representations of numbers can create an unacceptably large 
deviation between a computed binary representation and its theoretical numerical value due to 
successive rounding errors. This can be an intolerable situation when a higher degree of accuracy is 
20 required. Various methods can improve computation accuracy but they may require a significant 
increase in processing time and/or hardware. 

Agarwal et al in " New Scalar and Vector Elementary Functions for the IBM Svstem/370" 
IBM Journal of Research and Development: Vol.30 No.2 (March 1986), and Gal et al in "An 

CA9-2000-0040 I 



CA 0232S615 2000-11-10 



Accurate Elementary Mathematical Library for the IEEE Floating Point Standard " ACM 
Transactions on Mathematical Software Vol.17 No.l (March 1991) pp. 26 to 45, appear to disclose 
methods that include a lookup table having non-equidistantly spaced entries for computing values 
for elementary mathematical functions. 

5 Markstein in " Computation of Elementary Functions on the IBM RISC Svstem/6000 

Processor " IBM Journal of Research and Development Vol.34 No.l (January 1990) appears to 
disclose a method for computing values of a function g(x) near a given point x h 

Gustavson et al in " Elementary Function Subroutines" Department of Mathematical 
Sciences of IBM Watson Research Center (January 1985) appears to disclose a program 
10 implementation for a software mathematical library. 

Object of the Present Invention 

An object of the present invention is to obviate the disadvantages of the prior art. 

Summary of the Present Invention 

A first aspect of the present invention provides a machine-processing method for computing 

15 a property of a mathematically modelled physical system, the steps including: reading, via a 
machine processing unit, input data including a value for each identified ordered coefficient of a 
first polynomial p(x) representing the property, the polynomial p(x) being expressed as 
p(x) = • x J ) where j = 0 to «, a value of a quantity x, a value of a predetermined x h and a value 
of a predetermined pfc) correspondingly paired with said predetermined x f ; building, via the 

20 machine processing unit, a value of a second polynomial c(x) having ordered coefficients, the 
second polynomial c(x) being expressible as: c(x) = £(C* ■ **) where k = 0 to (w - 1 ) so that the first 
polynomial p(x) is expressible as: p(x) = p(Xi) + {x - jc,} • c(x), including the steps of: determining, 
via the machine processing unit, a value for each ordered coefficient of the second polynomial c(x) 
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by generating a plurality of machine processing unit signals to determine each of the ordered 
coefficients of the second polynomial c(x) from: C* = Z(iVi+/) -x^i) where j = 0 to (n - 1 - k) ; 
determining, via the machine processing unit, a value of the second polynomial c(x) by generating a 
plurality of machine processing unit signals to determine: c(x) = £(C* where k = 0 to (n - 1) ; 
5 and constructing, via the machine processing unit, a value of the first polynomial p(x) by generating 
a plurality of machine processing unit signals to determine:/?(x) =/?(x,) + {*-*/} • c(x); and 
outputting, via the machine-processing unit, the value of the first polynomial p(x) representing said 
property of the mathematically modelled physical system. 

Preferably the machine-implementable method is further adapted so that a difference 
io between x and x-, is sufficiently small to achieve a desired accuracy of a final computation of the 
machine representation of a numerical value of the first polynomial p(x). 

Preferably the machine-implementable method is further adapted so that the step of reading 
the input data includes reading, via the machine processing unit, the input data from a 
machine-readable medium. 

is Preferably the machine-implementable method is further adapted so that the ordered 

coefficients of the second polynomial c(x) are computed from a mathematical expression derivable 
from: C k = Z(iVi +J o • x*t) where j = 0 to (n - 1 - k). 

Preferably the machine-implementable method is further adapted so that the mathematical 
expression is a mathematical recurrence expression. 

20 Preferably the machine-implementable method is further adapted so that the mathematical 

recurrence expression is a forward mathematical recurrence expression. 
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Preferably the raachine-implementable method is further adapted so that the mathematical 
recurrence expression is a backward mathematical recurrence expression* 

Preferably the machine-implementable method is further adapted to implement the 
backward mathematical recurrence expression by including further steps for: equating, via the 
5 machine-processing unit, a value of a highest ordered coefficient of the second polynomial c(x) to a 
value of an identified highest ordered coefficient of the first polynomial p(x) by generating a 
plurality of machine processing unit signals to determine: CVi = P n \ and computing, via a machine 
processing unit, a value for each lower ordered coefficient of the second polynomial c(x) by 
generating a plurality of machine processing unit signals to determine: 
10 Ci-i = x f • Ck + Pk where /c=(n-l) to 1 . 

Preferably the machine-implementable method is further adapted so that the predetermined 
Xi is selected from a set of predetermined values of x f . 

Preferably the machine-implementable method is further adapted so that the predetermined 
Xi is a closest member of the set to the identified x. 

15 Preferably the machine-implementable method is further adapted so that the step of 

determining a value of the second polynomial c(x) is computed by using Horner's Rule. 

Preferably the machine-implementable method is further adapted for determining a value of 
a denominator polynomial q(x) having identified ordered denominator coefficients, the denominator 
polynomial q(x) being expressible as: q(x) = '*0 where ; = 0 to m, including further steps of: 
2 o adapting the input data to further include a value for each identified ordered denominator coefficient 
of the denominator polynomial q(x), a value of a predetermined qfa) correspondingly paired with 
the predetermined xr> and determining, via the machine processing unit, a value of another 
polynomial d(x) having ordered denominator coefficients, the another polynomial d(x) being 
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expressible as: d(x) = 2(Z)* where * = 0 to (m - 1 )so that the denominator polynomial q(x) is 
expressible as: q(x) = ?(*,) + {x-x f } *d(x) , and a value for the denominator polynomial is 
resolved. 

Preferably the machine-implementable method is further adapted so that the first 
5 polynomial p(x) is a numerator polynomial p(x), and p(x)-p(xj is computed, and p(x$ is not read. 

Preferably the machine-implementable method is further adapted for determining a value of 
a rational function r(x) being expressible as a quotient of the numerator polynomial p(x) and the 

/**) 

denominator polynomial q(x) expressed as r(x)= including further steps of: adapting the input 

data to further including a value of a predetermined r(x,) correspondingly paired with the 
10 predetermined x& constructing, via the machine processing unit, a value of the rational function r(x) 
by generating a plurality of machine processing unit signals to determine: 

q(xyq{Xj) p(xy-p{Xj) 
r(x) = r(x,).(l- ?(x) )+ 9(jc) . 

Preferably the machine-implementable method is further adapted so that the rational 
function is an approximation to a Bessel function. 

15 Preferably the machine-implementable method is further adapted so that the rational 

function is an approximation to an error function (ERF). 

Preferably the machine-implementable method is further adapted so that the rational 
function is an approximation to a complementary error function (ERFC). 

Preferably the machine-implementable method is further adapted so that the rational 
20 function is an approximation to a log gamma function (LGAMMA). 
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Preferably the machine-implementable method is further adapted so that all values are 
machine representations of some numerical value, the machine processing unit is a computer 
processing unit, and each machine representation is a binary representation operable with the 
computer processing unit, and the machine-readable medium is a computer-readable medium. 

A second aspect of the present invention provides a computer-readable program product 
having computer executable instructions for instructing a computer, the instructions tangibly 
embodying the machine-implementable method of the first aspect of the present invention. 

A third aspect of the present invention provides a computer-readable mathematical software 
routine library including computer executable instructions for instructing a computer, the 
instructions tangibly embodying the machine-implementable method of the first aspect of the 
present invention. 

Preferably, the computer-readable mathematical software routine library is further adapted 
so that the library is operatively associated with a software programming language. 

A fourth aspect of the present invention provides a machine for computing a property of a 
mathematically modelled physical system, the steps including: reading, via a machine processing 
unit, input data including a value for each identified ordered coefficient of a first polynomial p(x) 
representing the property, the polynomial p(x) being expressed as p(x) = £(Py • x J ) where j = 0 to «, 
a value of a quantity x, a value of a predetermined x h and a value of a predetermined p(x t ) 
correspondingly paired with the predetermined building, via the machine processing unit, a value 
of a second polynomial c(x) having ordered coefficients, the second polynomial c(x) being 
expressible as: c(x) = £(C* • x*)where k = 0 to (n - 1 )so that the first polynomial p(x) is expressible 
as: p(x) =p(xi)+ {x-Xi} 'C(x), including the steps of: determining, via the machine processing 
unit, a value for each ordered coefficient of the second polynomial c(x) by generating a plurality of 
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machine processing unit signals to determine each of the ordered coefficients of the second 
polynomial c(x) from: C* = • */) where j = 0 to (n - 1 - k)\ determining, via the machine 

processing unit, a value of the second polynomial c(x) by generating a plurality of machine 
processing unit signals to determine: c(x) = £(C* -x*) where k = 0 to (n - 1); constructing, via the 
machine processing unit, a value of the first polynomial p(x) by generating a plurality of machine 
processing unit signals to determine:p(x)=p(x/) + {x-x,} -c(x); and outputting, via the 
machine-processing unit, the value of the first polynomial p(x) representing the property of the 
mathematically modelled physical system. 

Preferably, the machine is further adapted so that a difference between x and x, is 
sufficiently small to achieve a desired accuracy of a final computation of the machine representation 
of a numerical value of the first polynomial p(x). 

Preferably, the machine is further adapted so that the means for reading the input data 
includes means for reading, via the machine processing unit, the input data from a machine-readable 
medium. 

Preferably, the machine is further adapted so that the ordered coefficients of the second 
polynomial c(x) are computed from a mathematical expression derivable from: 
C k = I,(P(m+j) • xft where j = 0 to (n - 1 - k) . 

Preferably, the machine is further adapted so that the mathematical expression is a 
mathematical recurrence expression. 

Preferably, the machine is further adapted so that the mathematical recurrence expression is 
a forward mathematical recurrence expression. 
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Preferably, the machine is further adapted so that the mathematical recurrence expression is 
a backward mathematical recurrence expression. 

Preferably, the machine is further adapted to implement the backward mathematical 
recurrence expression by further including: means for equating, via the machine processing unit, a 

5 value of a highest ordered coefficient of the second polynomial cfx) to a value of an identified 
highest ordered coefficient of the first polynomial pfx) by generating a plurality of machine 
processing unit signals to determine: Cn-\ = P n \ and means for computing, via the machine 
processing unit, a value for each lower ordered coefficient of the second polynomial cfx) by 
generating a plurality of machine processing unit signals to determine: 

10 Ck-i = xrCk+ Pk where k = (n- 1) to 1. 

Preferably, the machine is further adapted so that the predetermined x, is selected from a set 
of predetermined values of */. 

Preferably, the machine is further adapted so that the predetermined x t is a closest member 
of the set to the identified x. 

15 Preferably, the machine is further adapted for determining means for determining a value of 

the second polynomial c(x) is computed by using Horner's Rule. 

Preferably, the machine is further adapted for determining a value of a denominator 
polynomial q(x) having identified ordered denominator coefficients, the denominator polynomial 
qfx) being expressible as: q{x) = 2(0/ ■ x J ) where / = 0 to m, including further steps of: adapting the 
20 input data to further include a value for each identified ordered denominator coefficient of the 
denominator polynomial qfx), and a value of a predetermined qfxj correspondingly paired with the 
predetermined and determining, via the machine processing unit, a value of another polynomial 
dfx) having ordered denominator coefficients, the another polynomial dfx) being expressible as: 
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d(x) = • x*) where k = 0 to (m - 1 )so that the denominator polynomial q(x) is expressible as: 
q(x) = q(xi) + {* - Xi } • </(*), and a value for the denominator polynomial is resolved. 

Preferably, the machine is further adapted so that the first polynomial p(x) is a numerator 
polynomial p(x) 9 and p(x)-p(Xi) is computed, andpfxj is not read. 

5 Preferably, the machine is further adapted for determining a value of a rational function r(x) 

being expressible as a quotient of the numerator polynomial p(x) and the denominator polynomial 

P(x) 

q(x) expressed as r(x)= including further steps of: adapting the input data to further including 

a value of a predetermined r(x t ) correspondingly paired with the predetermined x,; and constructing, 
via the machine processing unit, a value of the rational function r(x) by generating a plurality of 

10 machine processing unit signals to determine: 

q(xyq(Xi) p(x)-p(Xj) 
r(x) = r{ Xi ) . (1 - q(pc) ) + q(x) . 

Preferably, the machine is further adapted so that the rational function is an approximation 
to a Bessel function. 

Preferably, the machine is further adapted so that the rational function is an approximation 
15 to an error function (ERF). 

Preferably, the machine is further adapted so that the rational function is an approximation 
to a complementary error function (ERFC). 

Preferably, the machine is further adapted so that the rational function is an approximation 
to a log gamma function (LGAMMA). 

20 Preferably, the machine is further adapted so that all values are machine representations of 

some numerical value, the machine processing unit is a computer processing unit, each machine 
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representation is a binary representation operable with the computer processing unit, and the 
machine-readable medium is a computer-readable medium. 

A fifth aspect of the present invention provides a machine having a computer-readable 
program product having computer executable instructions for instructing a computer to embody the 
5 machine of the fourth aspect of the present invention. 

A sixth aspect of the present invention provides a machine having a computer-readable 
mathematical software routine library including computer executable instructions for instructing a 
computer to embody the machine of the fourth aspect of the present invention. 

Preferably the computer-readable mathematical software routine library of the sixth aspect 
10 of the present invention is further adapted so that the library is operatively associated with a 
software programming language. 

Detailed Description of the Figures of the Present Invention 

To illustrate aspects of the present invention, the following figures are used, in which: 

Fig, 1 depicts a programming flow chart for computing a value of a polynomial p(x) 
15 in accordance with the present invention; 

Fig. 2 depicts a programming flow chart for computing a value of a rational function 
p(x)/q(x) in accordance with the present invention; and 

Fig. 3 depicts a programming flow chart for using a subroutine for computing a value 
of a Bessel function J 0 (x) in accordance with the present invention. 
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Detailed Description of a Preferred Embodiment of the Present Invention 

The present invention will be described with reference to an exemplary context of a 

computer processing method and apparatus for computing a binary representation of a numerical 

value of a polynomial p(x). 

5 For embodiments of the present invention, references to computing or calculating values or 

numbers will be understood to refer to programming instructions for instructing a computer to 
compute binary representations of numerical values of a mathematical expression or variable. It is 
understood that computing machinery can manipulate and transform binary representations of 
numerical values or numbers. 

10 A polynomial p(x) can be mathematically expressed as shown in the following equation: 

= £ Pjx j where./ = 0 tow in which p(x) = P 0 + Pt -x + P 2 -x 2 + ...+/>„ -x n Equ. 1 

Ordered coefficients of polynomial p(x) are Po 9 P n . Polynomial p(x) has /j+7 ordered 

coefficients and is said to have degree n. A numerical value of polynomial p(x) can be computed by 

reading each coefficient of polynomial p(x) and an identified x and using steps in Equation 1. 

n(n+\) 

15 Following Equation 1 directly would require a computer processor to perform 2 
multiplications and n additions. To significantly reduce the processing time, Horner's Rule can be 
applied to Equation 1, in which polynomial p(x) can be expressed in Equation la, as shown below: 

p(x) = P 0 +*(P 2 +*(...tfVi +* ■ Pn)-))) Equ. la 

By following Equation la, a program can perform n multiplications and n additions to 
20 compute a number for polynomial p(x) using less processing time than a program that follows 
Equation 1. 

A polynomial p(x) at x = x f can be expressed as follows: 
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pOci) = E Pjx*i wherey = 0tow 



Eqiu 2 



Therefore, a difference between Equation 1 and Equation 2 can be expressed as follows: 



p(x)-p(xi) = L Pj(* J -x?i) where; =0 ton 



Eqiu 3 



Expanding Equation 3 can provide the following expression: 



5 pC*)-/**,)- S {P/x-x,)-! [x^xi^lJwherey-Oton^-OtoO-l) 



Equ. 4 



For Equation 4, in the first sum ; = 0 to w, and in the second sum k = 0 to (j- 1 ). It is 
appreciated that Equation 4 can be expressed as follows: 

P(x) -p(x*) = (x- Xi ) . /Vi4,i -xft -x*} where 0 to(n- 1),; = 0 to(« - 1 -A:) 4a 

For Equation 4a, in the first sum k = 0 to(« - 1), and in the second sum / = 0 to(n - 1 - k). 
10 By following Equation 4a, a value for polynomial p(x) can be computed by reading coefficients of 
the polynomial (P J} P 2 , ... , P n \ and x, x„ and p(x$. Equation 4a can be expressed as a combination 
of Equation 5 and Equation 6, as follows: 

/t(x)-/7(x0= {x-x/}.c(x) ={x-x/}-E (C k -x k ) where A: = 0 to (n-1) Equ. 5 

C*-E/Vm*T*i where 7 = 0 to (n-\-k) Equ. 6 

is Equation 5 includes an expression for a new polynomial c(x) having coefficients, from C 0j 

Qwj, that are not known a priori and therefore need to be determined. Equation 6 can be used for 
determining values of the coefficients of new polynomial c(x) by involving the coefficients of 
polynomial p(x) and an identified x,. Optionally, other expressions for determining values for each 
coefficient of new polynomial c(x) can be derived from Equation 6. 

20 It is appreciated that Horner's Rule could be used to compute a value for polynomial c(x) 

after coefficients of polynomial c(x) have been computed by using Equation 6 or optionally from 
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another expression that could be derived from Equation 6. By using Horner's Rule, Equation 5 can 
be expressed as follows: 

p(x)-p(xi) = {x-xt} • {Co +x(Ci +x(C 2 +xt..(C„- 2 +*• C*.,)...)))} Equ. 5a 

Processing steps that follow Equation 5a begin within the innermost set of parentheses and 
proceed outwards to the outermost brackets. Calculated coefficients of new polynomial c(x) can 
subsequently be used in Equation 5a for computing a numerical value of polynomial p(x). 

Optionally, unknown coefficients of new polynomial c(x) can be computed by applying the 
principle of mathematical recurrence on Equation 6. For example, a forward mathematical 
recurrence can allow computation of coefficients of polynomial c(x) by beginning with the lowest 
ordered coefficient, Co, and proceeding to the highest ordered coefficient, C (n . ih Optionally, a 
backward mathematical recurrence can compute coefficients of polynomial c(x) by beginning with 
the highest ordered coefficient, C (n .jj and proceeding to the lowest ordered coefficient, C 0 . With each 
step of backwardly computing a coefficient of polynomial c(x) y a next outer most bracketed term of 
Equation 5a can be advantageously computed. 

Optionally, a backward recurrence for determining each coefficient of new polynomial c(x) 
can be realized by using Equation 7 and Equation 8. Each subsequent lower-ordered unknown 
coefficient of new polynomial c(x) can be computed by using Equation 8. 

Cn-\ = P n Equ. 7 

C k -\ = Xj-C k + Pk where (w- 1) to 1 Equ. 8 

To compute a numerical value of a polynomial p(x), a computer can read input data that 
includes a value for each identified ordered coefficient of the first polynomial p(x) 9 a value of a 
quantity x, a value of a predetermined x f , a value of a predetermined p(x$ that is correspondingly 
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paired with said predetermined x h An optional step can select a nearest x, that is closest to an 
identified x, in which the nearest x, is selected from a set of predetermined x h and a nearest 
predetermined p(x$ corresponds to the nearest X/. Optionally, coefficients of polynomial c(x) can be 
computed by following Equation 7 and Equation 8. A value for polynomial p(x) can then be 
5 computed by following Equation 5a. The operations of computing the coefficients of polynomial 
c(x) in Equation 8 can be interleaved with computing the nested parenthesis in Equation 5a so that 
each time a new Ct is computed by Equation 8 the next outer nested parenthesis of Equation 5a is 
computed next. 

Provided p(x) - p(x) does not have a multiple root at x = x„ the relative round-off error in 
10 computing p(x) - p(x) with Equation 5 can be kept small by keeping x sufficiently close to x,. This 
can be achieved by having available a number of sufficiently closely spaced x,'s covering the range 
of x to be handled, and using the nearest x f to x. To obtain a reasonable result that would be 
satisfactory for a given computational purpose, a difference between x and x, can be sufficiently 
small enough to achieve a desired accuracy of a computed value of polynomial p(x). 

is Optionally, accuracy can be further improved by choosing x, such that p(x) has extra 

accuracy beyond the precision of the floating-point number system of the computer. This can be 
achieved by selecting x, to be a floating-point number, by exhaustively searching from a starting 
desired value for x„ until an x, is found for which p(xj is very close to a floating-point number (i.e. 
closer than half the distance between adjacent floating-point numbers). Such a floating-point 

20 approximation to p(Xj) has an accuracy beyond the precision of the floating-point number system. 

Fig. 1 shows the preferred embodiment of the present invention, which is represented as a 
computer processing method embodied in a programming flow chart for computing a machine 
representation of a numerical value of a polynomial p(x). Optionally, coefficients ofp(x) can be 
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placed in a vector of length n. The method optionally includes a step of accessing a table having a 
set of predetermined values of x, and a corresponding set of predetermined values of p(x) each 
paired with a corresponding x,. A table can contain values for the following expressions: 

x/, p(Xi) for / = 1 to N. Equ. 8a 

In step 10, the method of the preferred embodiment begins. In step 15, each coefficient of 
polynomial p(x) and an identified x can be read by a computer processing unit. In step 20, the 
computer optionally selects a predetermined x t from a set of predetermined values of x,. Preferably, 
the predetermined x, is a closest member of the set to the identified x. In step 30, the initial 
conditions for a programming loop are set. Variable ksum is an ongoing summation that represents 
the sum in Equation Sa. Variable ksum is initially set equal to coefficient P n of polynomial p(x). 
Variable ck is a temporary value for any coefficients of polynomial c(x) 9 in which variable ck is 
computed to be a specific coefficient of polynomial c(x) that is determined in step 50 as will be 
shown below. Variable ck is initially set equal to P^ which is the highest ordered coefficient of 
polynomial p(x). Variable k is an index counter for a processing loop, in which variable k is initially 
set equal to (n-1), which is one less than the degree of polynomial p(x). In step 40, a query is 
performed to determine whether a value for variable k is equal to zero. If the numerical value for 
variable k is greater than zero, the method proceeds to step 50. If the numerical value for variable k 
is equal to zero, the method proceeds to step 60. In step 50, a value of variable ck is computed, 
which is coefficient C k by using the backward recurrence of Equation 8 that involves a higher 
ordered coefficient of polynomial p(x), which is P ky and a previously determined coefficient, Gw, of 
polynomial c(x). Next, an updated value for variable ksum is computed, which is an ongoing 
summation of the sum in Equation 5a, by determining a sum for the next outer bracketed term of 
Equation 5a. For remaining iterations of step 50, the summation represented in Equation 5a 
progresses from an innermost set of parentheses and proceeds outwardly for each iterative step in an 
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ongoing manner towards the outermost brackets. An updated value of variable k is computed so that 
the sum in Equation 5a can be resolved towards its outermost brackets. In step 60, the method uses 
a value for polynomial c(x), which is expressed in Fig. 1 as variable ksum after variable k equals 
zero in step 40. A value of polynomial p(x) is computed as a rearrangement of Equation 5a, which is 
5 expressed as p(x) = p(x) + (x-x$ *c(x). In step 70, the method for the preferred embodiment can stop 
or can be repeated from step 10. 

A second embodiment of the present invention provides a computer processing method for 
computing a machine representation of a numerical value of a rational function r(x). The rational 
function r(x) is expressed as a quotient of a numerator polynomial p(x) having a set of n 
10 predetermined coefficients from Po, Pi, P*> so that = £ Pj-x J where/ = 0 to n 9 and a 
denominator polynomial q(x) having m predetermined coefficients from Q 0 , Qu ...» Qm, so that 
#(*) = E where; = 0 torn. A rational function r(x) can be defined in the following 
expression: 



The second method can optionally include access to a table having a set of predetermined 
values of Xi and a corresponding set of predetermined values of q(x$ each paired with a 
corresponding x„ and a corresponding set of predetermined values of r(x$ each paired with a 
20 corresponding x h It is understood that each unique x, is paired with a correspondingly specific q(x)> 
and a correspondingly specific r(x), A spacing for table values of X/ can be set such that a resulting 
absolute value of (x-x<) is small enough to obtain a desired accuracy. The jc, spacing required for a 



r (*)=q(x) 



Equ. 8b 



15 



From Equation 8b, it follows that: 



q(x)-q(Xi) p(x)-p{xi) 
r(x) = r(x ; ).(l- q{x) )+ q(x) 



Equ. 9 
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given desired accuracy is typically estimated with a perturbation analysis (using Taylor series, for 
example), and verified by numerical testing. The closeness of x to *, needed to achieve a given 
accuracy in the computed r(x) depends on the particular rational function r(x). 

A table can contain values for the following expressions: 

5 x h r(x,), q{xd fori =1 to AT Equ. 10 

r(Xi) and q(xf) for i=ltoN Equ. 11 

Accuracy can be further improved by choosing the x's such that the following terms shown 
in Equation 11 can have extra accuracy beyond the precision of the floating-point number system of 
a computer. This can be achieved by selecting x f to be a floating-point number, by exhaustively 
10 searching from a starting desired value for x h until an *, is found for which so that both the 
quantities in Equation 11 are very close to floating-point numbers (closer than half the distance 
between adjacent floating-point numbers). Such floating-point approximations to r(x$ and q(x) have 
an accuracy beyond the precision of the floating-point number system. 

The second method computes values of a numerator polynomial p(x) and a denominator 
is polynomial q(x) by following Equation 5, as follows: 

P(x)-p(xi) = {jc-jc,} • c(x) Equ. 12a 

q(x)-q(xd = {*-*,} -d(x) Equ. 12 

Polynomial c(x) has (n-1) unknown coefficients and polynomial d(x) has (m-1) unknown 
coefficients. The second method can determine coefficients of c(x) and d(x) and then determine a 
20 value for rational function r(x). 

Fig. 2 shows the second embodiment of the present invention, which is represented as a 
computer processing method embodied in a programming flow chart for computing a computer 
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representation of a numerical value of a rational function r(x) = p(x)/q(x). Optionally, steps for 
accessing a table having a set of various predetermined values for x, a set of various predetermined 
values for q(x), and a set of various predetermined values for r(x$ can be included. 

In step 110, the second method begins. In step 115, a computer reads values of each 

5 coefficient of polynomial p(x) and each coefficient of polynomial q(x) y and an identified x. In step 
120, a predetermined x, from a set of predetermined values of x, can optionally be selected. 
Preferably, the predetermined x, is a closest member of the set to the identified x. These 
predetermined values can be optionally stored in a table or they can be supplied to the computer in 
some other means. Optionally, the table includes values for various xi, various q(xi) associated with 

10 the corresponding x;, and various r(x,) associated with the corresponding x h Optionally, specific 
values for the x/s, qfxfs, and rfxj's are chosen so that the q(xfs and r(x^s have extra accuracy 
beyond machine precision in the manner of the accurate table method that is disclosed by the prior 
art. It is appreciated that the computed result will be more accurate if the table has extra accuracy. 
Step 130, step 140, and step 150 are similar to step 30, step 40, and step 50 of Fig. 1, in that step 

is 130, step 140, and step 150 provide a value for variable hum to be used in step 160 as will be 
shown next. In step 160, a value for variable dp (i.e., delta p) is computed, which is a computed 
value for p(x)-p(x$ for a numerator polynomial p(x). Variable dp will be used in step 205 as will be 
shown later. Step 170, step 180, and step 190 are similar to step 130, step 140, and step 150 of Fig. 
2, in that step 170, step 180, and step 190 provide a value for variable ksum to be used in step 200 as 

20 will be shown next. In step 200 a value for variable dq (i.e., delta q), which is a computed value for 
q(x)-q(x$ for a denominator polynomial q(x) is computed. Variable dq will be used in step 202 and 
205 as shown later. In step 202 a numerical value for q(x) is computed by rearranging Equation 12. 
In step 205 a value for a rational function r(x) is computed by following Equation 9. In step 210 the 
second method can stop or can be restarted from step 110. 
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A third embodiment of the present invention provides a computer processing method for 
computing binary representations of a value for a Bessel function. A programming flow chart can be 
developed by adapting Fig. 2 for computing a value for a Bessel function Jo(x), which is a Bessel 
function of the first kind of order 0. The third method can be further adapted to compute a value for 
5 various Bessel functions, for example Jo, J/, Ko, and Yj y over a range of x (e.g., [0,8]) by first 
approximating the Bessel function by a piece- wise rational function. Bessel functions can typically 
be approximated by rational functions having numerator and denominator polynomials with the 
same degree (i.e. same number of coefficients). 

To compute a value of a Bessel function Jo(x) using the third method, the first task would be 
10 to determine a range that an identified x belongs in. The following is a list of steps for determining 
the proper range to which x belongs: if the identified x is a NAN (i.e., not a legitimate number, such 

as q"), or +/- INF (i.e., any non-zero number divided by zero), the method does not compute a value 
for the Bessel function; if the identified x is less than zero, the method makes the sign of the 
identified x positive; if the identified x is greater than or equal to zero and less than 3.72* 10* 9 , then 

15 the value of the Bessel function is 1, assuming the result is returned in IEEE double-precision 
floating point; if the identified x is greater than or equal to 3.72* 10 9 and less than 1, then the 
method uses a Taylor series to compute a value for the Bessel function; if the identified x is greater 
than or equal to 1 and less than 1.5, then the method calls a subroutine that incorporates the present 
invention for determining a value for the Bessel function. The subroutine will be explained in Fig. 

20 3, In each case where the subroutine is called, arguments P, Q, T 9 and n of the subroutine are read 
along with values of x and jc, associated with the particular range that x belongs with; if the 
identified x is greater than or equal to 1 .5 and less than 1 .9, then the method calls the subroutine for 
determining a value for the Bessel function; if the identified x is greater than or equal to 1.9 and less 
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than 2, then the method calls the subroutine for determining a value for the Bessel function; if the 
identified x is greater than or equal to 2 and less than 4, then the method performs steps that are 
beyond the scope of the present invention. The steps involve accurate evaluation of a rational 
approximation to J 0 (x)/(x-Zq) 9 where z 0 is the first zero of Jo\ if the identified x is greater than or 

5 equal to 4 and less than 4.5, then the method calls the subroutine for determining a value for the 
Bessel function; if the identified x is greater than or equal to 4.5 and less than 5, then the method 
calls the subroutine for determining a value for the Bessel function; if the identified x is greater than 
or equal to 5 and less than 6, then the method performs steps that are beyond the scope of the 
present invention. The steps involve accurate evaluation of a rational approximation to J 0 (x)/(x-z t ) 9 

10 where Zy is the first zero of J 0 \ if the identified x is greater than or equal to 6 and 6.4, then the 
method calls the subroutine for determining a value for the Bessel function; if the identified x is 
grater than 6.4 and less than 7, then the method calls the subroutine for determining a value for the 
Bessel function; if the identified x is greater than or equal to 7 and less than 7.5, then the method 
calls the subroutine for determining a value for the Bessel function; if the identified x is greater than 

is or equal to 7.5 and less than 8, then the method calls the subroutine for determining a value for the 
Bessel function; and if the identified x is greater than or equal to 8 then the method uses an 
asymptotic approximation for determining a value for the Bessel function, which is beyond the 
scope of the present invention. 

Fig. 3 shows the third embodiment of the present invention, which is represented as a 
20 computer processing method embodied in a programming flow chart for computing a machine 
representation of a numerical value of a Bessel function J 0 (x). 

In step 360 a subroutine, which can be called eval_rat, of the flow chart of the third 
embodiment begins. The subroutine inputs an argument x, at which the value of a rational function 
is desired to be computed, and P and Q, which are vectors containing the coefficients of the 
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numerator and denominator polynomials p(x) and q(x), respectively, of the rational function, and 7*, 
which is a table containing a set of predetermined x/s and corresponding r(x$ 's and q(xfs, and w, 
which is the degree of p(x) and q(x). In step 362 a computer processing unit reads values for each 
coefficient of numerator polynomial p(x) and denominator polynomial q(x)> in which both 
5 polynomials have the same number or coefficients, and also reads an identified x. In step 364, to 
achieve improved performance, a nearest X/ and other required data are optionally read from a table. 
In step 366 various variables are initialized for subsequent use in an instruction loop. One software 
loop is used because the numerator and the denominator polynomials have the same degree (that is, 
share the same number of coefficients). This loop is similar to a combination of the two loops 
10 involving steps 130 to ISO and 170 to 190 in Fig, 2. In step 368, a test determines whether the loop 
has been completed. In step 370, variables are computed in a manner similar to that combining 
steps 150 and 190 in Fig. 1. In step 372, variables are computed in a manner similar to that 
combining steps 160 and 200 to 205 shown in Fig. 1. Step 374 marks the end of the subroutine of 
the flow chart of the third embodiment. 

15 For non-elementary special functions, for example the error function (ERF), the 

complementary error function (ERFC), or the Bessel functions,range-reduction properties, such as 
those available for the elementary functions, are not known. To directly apply prior art table-driven 
techniques, for example those used for the elementary functions, for these non-elementary special 
functions, would require a computer processing unit to read coefficients of each of a large number 

20 of polynomials, one associated with each table point. This can make the table unacceptably large. 

Computing steps that embody the methods of the present invention can be used to 
determine values for polynomials or rational functions approximating the non-elementary special 
function. Optionally, coefficients of a polynomial or rational function associated with each 
individual table point, x i9 do not have to be stored; instead, these coefficients can be computed using 
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simple formulas from only one (for polynomials) or two (for rational functions) stored function 
values per table point, resulting in tables of manageable size. In this sense, the present invention 
provides a benefit similar to that of range reduction formulas, for functions for which range 
reduction formulas do not exist. 

Mathematical software libraries, including those associated with a software programming 
language or operating system, can be adapted to instruct a general purpose digital computer to 
embody the present invention. It is appreciated that a software does not necessarily have to belong 
to a programming language or operating system. There are also stand-alone libraries of 
mathematical routines, for example, the IBM™ Engineering and Scientific Subroutine Library, 
ESSL™. 

There are many examples that illustrate the practical application of the present invention for 
mathematically modelling a property of a physical system. C M. Close in " The Analysis of Linear 
Circuits ". Harcourt, Brace & World (1966), discloses in Chapters 6 and 11 the use of polynomials 
for determining frequency response characteristics of electrical systems. The same concepts can be 
applied to non-electrical systems, such as mechanical, hydraulic, acoustical, and thermal systems. 
R.C. Weyrick in " Fundamentals of Automatic Control ". McGraw-Hill Book Company (1975), 
discloses in Chapters 2 and 6 the use of polynomials for mathematically modelling physical systems 
for predicting a behavior of a physical system. D. Halliday and R. Resnick in " Physics: Parts 1 and 
2" John Wiley & Sons (1978) disclose in Chapters 3 and 4 the use of polynomials for 
mathematically modelling and predicting spatial co-ordinates of impact of a projectile being 
released from an aeroplane. 

The present invention can be implemented in a computer-readable computer memory 
program having a medium for storing and carrying the program to a general purpose computer for 
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instructing a computer processing unit to implement the embodiments of the present invention. A 
computer-memory product can be a floppy disk or compact disk that can be adapted to carry 
software that tangibly embodies the present invention. 

A computer processing method that tangibly embodies the present invention can be used for 
computing a number for a mathematical model and for comparing the computed number against an 
observed behaviour of a physical phenomenon. 

The concepts of the present invention can be further extended to a variety of other 
applications that are clearly within the scope of this invention. 

Having thus described the present invention with respect to a preferred embodiment as 
implemented, it will be apparent to those skilled in the art that many modifications and 
enhancements are possible to the present invention without departing from the basic concepts as 
described in the preferred embodiment of the present invention. Therefore, what is intended to be 
protected by way of letters patent should be limited only by the scope of the following claims. 
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The embodiments of the invention in which an exclusive property or privilege is claimed are 
defined as follows: 

1. A machine-processing method for computing a property of a mathematically modelled 
5 physical system, the steps comprising: 

a) reading, via a machine processing unit, input data including a value for each identified 
ordered coefficient of a first polynomial p(x) representing said property, said polynomial p(x) being 
expressed as p(x) = X(P, • x J ) where j = 0 to n 9 a value of a quantity a value of a predetermined 
x i9 and a value of a predetermined p(x$ correspondingly paired with said predetermined xr 9 

10 b) building, via said machine processing unit, a value of a second polynomial c(x) having 

ordered coefficients, said second polynomial c(x) being expressible as: 
c(x) = • x*) where k=0 to (h-I)so that said first polynomial p(x) is expressible as: 
pipe) = p(pd) + {x -x f ) • c(x), comprising the steps of: 

i) determining, via said machine processing unit, a value for each ordered coefficient 
15 of said second polynomial c(x) by generating a plurality of machine processing unit signals to 

determine each said ordered coefficient of said second polynomial c(x) from: 
C* = E(/W,v*i) where ; = 0 to («- 1 -*) ; 

ii) determining, via said machine processing unit, a value of said second polynomial 
c(x) by generating a plurality of machine processing unit signals to determine: 

20 c(x) = >x k ) where k= 0 to (n- 1) ; 
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c) constructing, via said machine processing unit, a value of said first polynomial p(x) by 
generating a plurality of machine processing unit signals to determine :p(x) =p(xi) + -c(x); 
and 

d) outputting, via said machine-processing unit, said value of the first polynomial p(x) 
5 representing said property of the mathematically modelled physical system. 

2. The machine-implementable method of claim 1, wherein a difference between x and x* is 
sufficiently small to achieve a desired accuracy of a final computation of said machine 
representation of a numerical value of said first polynomial p(x). 

10 

3. The machine-implementable method of claim 2 wherein the step of reading said input data 
comprises reading, via said machine processing unit, said input data from a machine-readable 
medium. 

15 4. The machine-implementable method of claim 3 wherein said ordered coefficients of said 
second polynomial c(x) are computed from a mathematical expression derivable from: 
C* = E(JVi*r*i) where 7 = 0 to (n-l-k). 

5. The machine-implementable method of claim 4 wherein said mathematical expression is a 
2 o mathematical recurrence expression. 
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6, The machine-implementable method of claim 5 wherein said mathematical recurrence 
expression is a forward mathematical recurrence expression. 

7, The machine-implementable method of claim 5 wherein said mathematical recurrence 
5 expression is a backward mathematical recurrence expression. 

8. The machine-implementable method of claim 7 further adapted to implement said backward 
mathematical recurrence expression by comprising further steps for: 

e) equating, via said machine-processing unit, a value of a highest ordered coefficient of said 
10 second polynomial c(x) to a value of an identified highest ordered coefficient of said first 

polynomial p(x) by generating a plurality of machine processing unit signals to determine: 
Cn-\ = Pn\ and 

f) computing, via a machine processing unit, a value for each lower ordered coefficient of 
said second polynomial c(x) by generating a plurality of machine processing unit signals to 

is determine: C*-i = x, • C* + Pk where k = (n- 1) to 1. 

9. The machine-implementable method of claim 8 wherein said predetermined x f is selected 
from a set of predetermined values of x,% 

20 10. The machine-implementable method of claim 9 wherein said predetermined Xj is a closest 
member of said set to said identified x. 
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11. The machine-implementable method of claim 10 wherein said step of determining a value of 
said second polynomial cfx) is computed by using Horner's Rule. 

5 12. The machine-implementable method of claim 1 1 for determining a value of a denominator 
polynomial q(x) having identified ordered denominator coefficients, said denominator polynomial 
q(x) being expressible as: q(x) = »*0 where j = 0 to m, comprising further steps of: 

g) adapting said input data to further include a value for each identified ordered denominator 
coefficient of said denominator polynomial qfx), a value of a predetermined qfx) correspondingly 

10 paired with said predetermined x»\ and 

h) determining, via said machine processing unit, a value of another polynomial d(x) having 
ordered denominator coefficients, said another polynomial dfx) being expressible as: 
d(x) = 2(0* • x*) where k = 0 to (m - 1 )so that said denominator polynomial qfx) is expressible as: 
q(x) = q(xj) + {x - Xj } - dfx) , and a value for the said denominator polynomial is resolved. 

15 

13. The machine-implementable method of claim 12 wherein the first polynomial pfx) is a 
numerator polynomial pfx), and pfx)-pfx) is computed, and pfx) is not read. 



14. The machine-implementable method of claim 13 for determining a value of a rational 
20 function rfx) being expressible as a quotient of said numerator polynomial pfx) and said 

P(x) 

denominator polynomial qfx) expressed as rfx)- , comprising further steps of: 
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i) adapting said input data to further including a value of a predetermined r(x/) 
correspondingly paired with said predetermined x,; and 

j) constructing, via said machine processing unit, a value of said rational function r(x) by 
generating a plurality of machine processing unit signals to determine: 

q(x)-g(Xi) p(xyp(Xj) 
5 r(x) = r(x,).(l- ) + q{pc ) • 

15. The machine-implementable method of claim 14 wherein said rational function r(x) is an 
approximation to a Bessel function. 

10 16. The machine-implementable method of claim 14 wherein said rational function r(x) is an 
approximation to an error function (ERF). 

17. The machine-implementable method of claim 14 wherein said rational function r(x) is an 
approximation to a complementary error function (ERFC). 

15 

18. The machine-implementable method of claim 14 wherein said rational function r(x) is an 
approximation to a log gamma function (LGAMMA). 

19. The machine-implementable method of claim 11 or 14 wherein all values are machine 
20 representations of some numerical value, said machine processing unit is a computer processing 
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unit, each machine representation is a binary representation operable with said computer processing 
unit, and machine-readable medium is a computer-readable medium. 

20. A computer-readable program product having computer executable instructions for 
instructing a computer, said instructions tangibly embodying the machine-implementable method of 
claim 19. 

21. A computer-readable mathematical software routine library including computer executable 
instructions for instructing a computer, said instructions tangibly embodying the 
machine-implementable method of claim 19. 

22. The computer-readable mathematical software routine library of claim 21 wherein said 
library is operatively associated with a software programming language. 
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23. A machine for computing a property of a mathematically modelled physical system, the 
steps comprising: 

a) reading, via a machine processing unit, input data including a value for each identified 
ordered coefficient of a first polynomial p(x) representing said property, said polynomial p(x) being 
expressed as p(x) = • x J ) where j = 0 to n> a value of a quantity x, a value of a predetermined 
*i, and a value of a predetermined p(Xj) correspondingly paired with said predetermined xr, 

b) building, via said machine processing unit, a value of a second polynomial c(x) having 
ordered coefficients, said second polynomial c(x) being expressible as: 
c(x) = 2 t (Ck'X k ) where A = 0 to (w-l)so that said first polynomial p(x) is expressible as: 
p(x) = p(xi) + {x- Xi} • c(x), comprising the steps of: 

i) determining, via said machine processing unit, a value for each ordered coefficient 
of said second polynomial c(x) by generating a plurality of machine processing unit signals to 
determine each said ordered coefficient of said second polynomial c(x) from: 
C* = Z(Pot + t + W/) where ; = 0to(«-l-*); 

ii) determining, via said machine processing unit, a value of said second polynomial 
c(x) by generating a plurality of machine processing unit signals to determine: 
c(x) = £(C* • x k ) where A: = 0to (n- 1); 

c) constructing, via said machine processing unit, a value of said first polynomial p(x) by 
generating a plurality of machine processing unit signals to determine:/?(x) =p(x,) + {x-x,} -c(x)\ 
and 

d) outputting, via said machine-processing unit, said value of the first polynomial p(x) 
representing said property of the mathematically modelled physical system. 
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24. The machine of claim 23 wherein a difference between x and x* is sufficiently small to 
achieve a desired accuracy of a final computation of said machine representation of a numerical 
value of said first polynomial p(x). 

5 

25. The machine of claim 24 wherein said means for reading said input data comprises means 
for reading, via said machine processing unit, said input data from a machine-readable medium. 

26. The machine of claim 25 wherein said ordered coefficients of said second polynomial c(x) 
10 are computed from a mathematical expression derivable from: 

C k = J! / (P(M+j)-x f i ) where ; = 0 to (n-\-k). 

27. The machine of claim 28 wherein said mathematical expression is a mathematical recurrence 
expression. 

15 

28. The machine of claim 27 wherein said mathematical recurrence expression is a forward 
mathematical recurrence expression. 

29. The machine of claim 27 wherein said mathematical recurrence expression is a backward 
20 mathematical recurrence expression. 
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30. The machine of claim 29 further adapted to implement said backward mathematical 
recurrence expression by further comprising: 

e) means for equating, via said machine processing unit, a value of a highest ordered 
coefficient of said second polynomial c(x) to a value of an identified highest ordered coefficient of 
said first polynomial p(x) by generating a plurality of machine processing unit signals to determine: 
Cn-\ = P n \ and 

f) means for computing, via said machine processing unit, a value for each lower ordered 
coefficient of said second polynomial c(x) by generating a plurality of machine processing unit 
signals to determine: Ck-\ = *, • C* + P* where k = (n - 1) to 1 . 

31. The machine of claim 30 wherein said predetermined x- t is selected from a set of 
predetermined values of*,. 

32. The machine of claim 30 wherein said predetermined x, is a closest member of said set to 
said identified x. 

33. The machine of claim 32 wherein the determining means for determining a value of said 
second polynomial c(x) is computed by using Horner's Rule. 
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34. The machine of claim 33 for determining a value of a denominator polynomial q(x) having 
identified ordered denominator coefficients, said denominator polynomial q(x) being expressible as: 
qix) = 2(0/ * *0 where j = 0 to m, comprising further steps of: 

g) adapting said input data to further include a value for each identified ordered denominator 
5 coefficient of said denominator polynomial q(x) 9 and a value of a predetermined qfxj 

correspondingly paired with said predetermined and 

h) determining, via said machine processing unit, a value of another polynomial d(x) having 
ordered denominator coefficients, said another polynomial d(x) being expressible as: 
d(x) = E(£>* • x k ) where k= 0 to (m - 1 )so that said denominator polynomial q(x) is expressible as: 

io q(x) = q(xj) + {x - Xi } ■ d(x) , and a value for the said denominator polynomial is resolved. 

35. The machine of claim 34 wherein the first polynomial p(x) is a numerator polynomial p(x), 
and p(x)-p(xi) is computed, and p(X() is not read. 

is 36. The machine of claim 35 for determining a value of a rational function r(x) being 
expressible as a quotient of said numerator polynomial p(x) and said denominator polynomial q(x) 
P(x) 

expressed as r(x)= , comprising further steps of: 

i) adapting said input data to further including a value of a predetermined r(x,) 
correspondingly paired with said predetermined jc,; and 

20 j) constructing, via said machine processing unit, a value of said rational function r(x) by 

generating a plurality of machine processing unit signals to determine: 
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g(xyq(Xj) p(xy-p(Xj) 
r(x) = r(x,)-(l- q(x) ) + q(x) 

37. The machine of claim 36 wherein said rational function is an approximation to a Bessel 
function. 

5 

38. The machine of claim 36 wherein said rational function is an approximation to an error 
function (ERF). 

39. The machine of claim 36 wherein said rational function is an approximation to a 
10 complementary error function (ERFC). 

40. The machine of claim 36 wherein said rational function is an approximation to a log gamma 
function (LGAMMA). 

is 41, The machine of claim 33 or 16 wherein all values are machine representations of some 
numerical value, said machine processing unit is a computer processing unit, each machine 
representation is a binary representation operable with said computer processing unit, and said 
machine-readable medium is a computer-readable medium. 
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42. A machine having a computer-readable program product having computer executable 
instructions for instructing a computer to embody the machine of claim 41 . 

43. A machine having a computer-readable mathematical software routine library including 
computer executable instructions for instructing a computer to embody the machine of claim 41. 

44. A machine having the computer-readable mathematical software routine library of claim 43 
wherein said library is operatively associated with a software programming language. 
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Qo. Qi Q m , and an 




identified x. 



t 



Read a nearest x, from a 
predetermined set of x,'s. 



120 



200. 



202. 



205 



Set: ksum = P n 
Set: ck = P n 
Set k = n-l 



130 




160 


Compute: 




dp = (x-x,) • ksum 




t 


170 — 


Set: ksum = Q m 






■ Set:ck=Q m 






Set: k = m-1 





Compute: 
ck = P k +(x, *ck) 
ksum = ck + (x * ksum) 
k = k-1 




190 



Compute: 
dq = (x - x ( ) * ksum 

x 



Compute: 
qx°<yQ + dq 



Compute: 
ck = Q k +(x,*ck) 
ksum = ck + (x* ksum) 
k = k-1 



Return the Result: 
oc = fx*P| - do] +dE 
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Figure 3 



372 




START 
evaUat (x, P, Q, T. n) 




I 



Read 

Po, Pt, .... Pn, 

Qo, Qi Q„,and 

an identified x. 



Select and present a nearest 
Xi from a set of predetermined 
x/s such that the nearest x, is 
closest to the identified x. 



Set ksump = P n , ksumq = Q n 
Set: ckp = P n ,ckq = Q n 
Set: k = n-1 
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364 



,366 



NO 



370 





Compute: 




dp = (x - x,) * ksump 




dq = (x - Xj) * ksumq 




qx = x,+ dq 




Return the Result: 




nc = rxi+(dp-rx,*dq)/qx 
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Compute: 
ckp= P t +(x*ckp) 
ksump = ckp + (x * ksump) 

ckq= Q k +(x»*ckq) 
ksumq = ckq + (x * ksumq) 
k = k-1 



